Graphene¶
Graphene consists of carbon atoms in a 2d hexagonal configuration. Like WSe\(_2\), graphene is a popular material for twisting and can be considered the cradle of twistronics since Macdonald’s seminal publication. [BM11] Due to its simplicity and similarity to our research interest, it is a good example to showcase the application of Quantum Espresso (QE) [GAB+17] and Wannier90 (W90) [PVA+20] to a 2d material.
Basic theory¶
A possible choice of lattice vectors for graphene is:
where length is in units of the lattice spacing \(a\) between two carbon atoms, which is approximately 1.42 Å. [CDAnjouG+11]
Fig. 13 Lattice structure of graphene.¶
In the case of sublattice symmetry the Hamiltonian writes:
with:
where energy is measured in units of the hopping parameter \(t\approx 2.7\) eV between two carbon atoms [CDAnjouG+11].
Fig. 14 Tight-binding bandstructure of graphene from (22). At the K point there is a Dirac point with linear dispersion.¶
In Fig. 14 we plot the band structure for graphene for this tight-binding Hamiltonian. At the K point the dispersionbecomes linear and forms a Dirac point.
Quantum Espresso¶
Now we wish to recreate this bandstructure and Hamiltonian with a dft calculation. The first step here is a self-consistent field calculation which fixes the electron density for future band structure calculations.
SCF/NSCF/Bands input¶
scf &CONTROL input
&CONTROL
prefix = 'graphene'
outdir = './out'
verbosity = 'high'
pseudo_dir = '../pseudo/'
calculation = 'nscf'
/
Description
The &CONTROL card contains parameters controlling IO
QE finds the required pseudopotentials in the
pseudo_dirand stores them and all datafiles in theoutdir. Theprefixfunctions as a label for these files so that later calculations can find these specific files and reuse them.verbositysets the amount of information to be written to the output file.calculationtype, in this case scf.
scf &SYSTEM input
&SYSTEM
assume_isolated = '2D'
ibrav = 4
nat = 2
ntyp = 1
occupations = 'smearing'
degauss = 0.015
ecutwfc = 45
nbnd = 30
a = 2.45951214
c = 32
/
Description
The &SYSTEM card sets parameters controlling lattice properties
QE contains a number of preset lattices, with
ibrav4 we choose a triangular lattice with lattice vectors:\(v_1 = a\begin{pmatrix}1&0&0\end{pmatrix}\)
\(v_2 = a\begin{pmatrix}-\frac12&\frac12\sqrt3&0\end{pmatrix}\)
\(v_3=a\begin{pmatrix}0&0&\frac ca\end{pmatrix}\)
Here
aandccontroll the size of the lattice. For 2d lattices, like graphene, we useassume_isolatedand a large value forcto make sure there is no periodicity in the z direction. Finallynatandntypeinform QE of the total number and the number of unique atoms in our unit cell.In the case of metals DFT can be unstable while integrating the Fermi surface. We prevent problems by setting
occupationstosmearingwhich smears the bands with a value set bydegauss.ecutwfcis a cut-off for the plane wave expansion of the wavefunctions. We can find appropriate cut-offs for each atom type at the materialscloud. [PMC+18]nbnbis the number of bands we wish QE to calculate. For a complete picture of the conduction bands it is important to set this parameter manually.
scf &ELECTRONS input
&ELECTRONS
/
Description
The &ELECTRONS card must be included for a scf calculation even though we use only the default parameters.
scf ATOMIC inputs
ATOMIC_SPECIES
C 12.0107 'C.pbe-n-kjpaw_psl.1.0.0.UPF'
ATOMIC_POSITIONS angstrom
C 0.0000000 0.0000000 0.0000000
C 0.0000000 1.4200000 0.0000000
Description
The
ATOMIC_SPECIEScard is a list of all atom types with their weight and pseudopotential files, non-relativistic pseudopotentials can also be found on the materialscloud. [PMC+18]The
ATOMIC_POSITIONScard sets the unit of length and the positions of the atoms in the unit cell. We use angstrom as our unit of choice.
scf K_POINTS input
K_POINTS automatic
9 9 1 1 1 1
Description
The K_POINTS card sets the k points at which QE calculates the wavefunctions and energies. The automatic option with input \(n_1, n_2,n_3,s_1, s_2, s_3\) generates a Monkhorst-Pack grid of dimension \(n_1\times n_2\times n_3\) with respective offsets \(s_i=0\vee1\), corresponding to no offset or half a gridpoint offset.
The bands and nscf files are identical save for the calculation and K_POINTS cards:
nscf K_POINTS input
K points list
K_POINTS crystal
64
0.0000000 0.0000000 0.0000000 1
0.0000000 0.1250000 0.0000000 1
0.0000000 0.2500000 0.0000000 1
0.0000000 0.3750000 0.0000000 1
0.0000000 0.5000000 0.0000000 1
0.0000000 0.6250000 0.0000000 1
0.0000000 0.7500000 0.0000000 1
0.0000000 0.8750000 0.0000000 1
0.1250000 0.0000000 0.0000000 1
0.1250000 0.1250000 0.0000000 1
0.1250000 0.2500000 0.0000000 1
0.1250000 0.3750000 0.0000000 1
0.1250000 0.5000000 0.0000000 1
0.1250000 0.6250000 0.0000000 1
0.1250000 0.7500000 0.0000000 1
0.1250000 0.8750000 0.0000000 1
0.2500000 0.0000000 0.0000000 1
0.2500000 0.1250000 0.0000000 1
0.2500000 0.2500000 0.0000000 1
0.2500000 0.3750000 0.0000000 1
0.2500000 0.5000000 0.0000000 1
0.2500000 0.6250000 0.0000000 1
0.2500000 0.7500000 0.0000000 1
0.2500000 0.8750000 0.0000000 1
0.3750000 0.0000000 0.0000000 1
0.3750000 0.1250000 0.0000000 1
0.3750000 0.2500000 0.0000000 1
0.3750000 0.3750000 0.0000000 1
0.3750000 0.5000000 0.0000000 1
0.3750000 0.6250000 0.0000000 1
0.3750000 0.7500000 0.0000000 1
0.3750000 0.8750000 0.0000000 1
0.5000000 0.0000000 0.0000000 1
0.5000000 0.1250000 0.0000000 1
0.5000000 0.2500000 0.0000000 1
0.5000000 0.3750000 0.0000000 1
0.5000000 0.5000000 0.0000000 1
0.5000000 0.6250000 0.0000000 1
0.5000000 0.7500000 0.0000000 1
0.5000000 0.8750000 0.0000000 1
0.6250000 0.0000000 0.0000000 1
0.6250000 0.1250000 0.0000000 1
0.6250000 0.2500000 0.0000000 1
0.6250000 0.3750000 0.0000000 1
0.6250000 0.5000000 0.0000000 1
0.6250000 0.6250000 0.0000000 1
0.6250000 0.7500000 0.0000000 1
0.6250000 0.8750000 0.0000000 1
0.7500000 0.0000000 0.0000000 1
0.7500000 0.1250000 0.0000000 1
0.7500000 0.2500000 0.0000000 1
0.7500000 0.3750000 0.0000000 1
0.7500000 0.5000000 0.0000000 1
0.7500000 0.6250000 0.0000000 1
0.7500000 0.7500000 0.0000000 1
0.7500000 0.8750000 0.0000000 1
0.8750000 0.0000000 0.0000000 1
0.8750000 0.1250000 0.0000000 1
0.8750000 0.2500000 0.0000000 1
0.8750000 0.3750000 0.0000000 1
0.8750000 0.5000000 0.0000000 1
0.8750000 0.6250000 0.0000000 1
0.8750000 0.7500000 0.0000000 1
0.8750000 0.8750000 0.0000000 1
Description
The K_points card is list of points along a path in k space for bands and a grid for nscf, preceded by the number of points, in this case 64. Here we replace automatic with crystal which reads the k points in crystal coordinates (relative to the reciprocal lattice vectors). The fourth number is the relative weight of the k point.
Command line interface¶
The command line interface for QE is quite simple:
pw.x < prefix.scf.in > prefix.scf.out
where prefix is graphene in our case.
We can parallelize QE with the mpi library:
mpirun -n #n pw.x < prefix.scf.in > prefix.scf.out
here the -n flag can differ per distribution (sometimes it is -np). #n is the number of processes at our disposal. We can then choose how to divide these processes with additional flags after the pw.x executable. The main flag we used is the -nk flag which divides the calculation in #nk pools which each calculates a fraction of the k points, allowing for perfect parallelization. For higher #nk the RAM usage also increases so we can not automatically just set #nk to the amount of k points.
Bands output¶
The main output we directly care about is the bands output which leads a summary of the input arguments and the symmetries of the lattice. Since we define all K points explicitly these will not be used by QE. When we choose an automatic grid QE significantly decreases the points it needs to evaluate. Near the end of the file all k points and energies of the bands are printed.
bands output
k = 0.0000 0.0000 0.0000 ( 5649 PWs) bands (ev):
-23.8181 -11.9014 -7.2766 -7.2766 -1.3477 -0.9995 -0.7246 -0.5157
-0.0219 0.2830 1.0299 1.4288 2.4049 2.9177 4.0662 4.1555
4.1555 4.7234 6.0004 6.8461 7.2901 8.1773 8.4570 9.2604
10.5975 11.9548 13.2678 14.9152 16.2216 18.1415
Description
The first k point evaluated by QE and all corresponding energies in eV.
Fig. 15 QE bandstructure of graphene with the bands from the TB Hamiltonian included, translated so that the K points overlap. The TB bands match the valence bands exceptionally well.¶
In Fig. 15 we plot the QE bands together with the TB bands from (22). The model works quite well for the valence part of the QE bands. Ostensibly we require NNN hoppings and beyond to capture the conduction band component more accurately.
Projwfc input¶
Before we start the wannierization process we need to know which bands we want to wannierize and which atomic orbital contribute to these bands. QE offers a program called projwfc (project wavefunctions) which does exactly this.
projwfc input
&PROJWFC
prefix = 'graphene'
outdir = '../out'
DeltaE = 0.05
/
Description
The PROJWFC card is the only card for projwfc source files. The prefix and outdir should match the project. DeltaE is the step size for the dos calculation of the projected orbitals, which we do not need but must be set for the program to work.
Projwfc output¶
PROJWFC output: states
state # 1: atom 1 (C ), wfc 1 (l=0 m= 1)
state # 2: atom 1 (C ), wfc 2 (l=1 m= 1)
state # 3: atom 1 (C ), wfc 2 (l=1 m= 2)
state # 4: atom 1 (C ), wfc 2 (l=1 m= 3)
state # 5: atom 2 (C ), wfc 1 (l=0 m= 1)
state # 6: atom 2 (C ), wfc 2 (l=1 m= 1)
state # 7: atom 2 (C ), wfc 2 (l=1 m= 2)
state # 8: atom 2 (C ), wfc 2 (l=1 m= 3)
Description
The projwfc output file contains a list of the projected atomic orbitals. For graphene we have two sites: atom 1 and atom 2, two atomic orbital types: s-orbitals with azimuthal quantum number l=0 (wfc 1) and p-orbitals with azimuthal quantum number l=1 (wfc 2). For l=1 there are three different possibilities for the magnetic quantum number \(m_l=-1, 0, 1\). For a spinless bandstructure this quantum number is reformulated to m=1 for the p\(_z\) orbital, m=2 for the p\(_x\) orbital and m=3 for the p\(_y\) orbital. In total we have 8 orbitals.
The second part of the projwfc output file list the contribution of each state to each band at each k points.
PROJWFC output: states
k = 0.0000000000 0.0000000000 0.0000000000
==== e( 1) = -23.81812 eV ====
psi = 0.492*[# 1]+0.492*[# 5]+
|psi|^2 = 0.984
==== e( 2) = -11.90144 eV ====
psi = 0.473*[# 2]+0.473*[# 6]+
|psi|^2 = 0.946
==== e( 3) = -7.27665 eV ====
psi = 0.249*[# 3]+0.249*[# 4]+0.249*[# 7]+0.249*[# 8]+
|psi|^2 = 0.997
==== e( 4) = -7.27665 eV ====
psi = 0.249*[# 3]+0.249*[# 4]+0.249*[# 7]+0.249*[# 8]+
|psi|^2 = 0.997
==== e( 5) = -1.34766 eV ====
psi = 0.002*[# 1]+0.002*[# 5]+
|psi|^2 = 0.005
==== e( 6) = -0.99950 eV ====
psi = 0.004*[# 2]+0.004*[# 6]+
|psi|^2 = 0.007
Description
The contribution of the eight states to the first 6 bands at the first k point. If we write the wavefunction as \(\Psi = \sum\alpha_i\psi_i\) then the number before each state corresponds to \(|\alpha_i|^2\). We see that for bands 5 and 6 |psi|^2 is not even close to 1. This tells us that this band can not be described with these 8 states. The leftover plus at the end of the line is just a for loop artifact which can be ignored.
In tandem with the projected orbitals we wrote a small script to untangle the bands so we can separate the future Wannier bands from the other conduction bands.
Fig. 16 QE bandstructure of graphene with possible Wannier bands in solid black, left over bands are dotted.¶
Wannier90¶
Next up is the wannierization of the two bands forming the Dirac cone. FromFig. 16 we learn that only the p\(_z\) orbitals contribute and that these bands extend far into the valence bands. To untangle these bands we need to make sure that takes all these bands into account. That is why we have set nbnb quite high for the nscf calculation so that QE produces the whole band. It is essential that the nscf calculation is run after a bands calculation using the same outdir for both as they overwrite each other’s wavefunctions.
W90 input¶
W90 general input
guiding_centres = true
write_hr = true
wannier_plot = true
bands_plot = true
num_bands = 30
num_wann = 2
num_iter = 2000
wannier_plot_supercell = 3, 3, 1
Description
guiding_centrestells W90 to use the projection (orbital) centres as the guiding centres for the Wannier orbitalswrite_hrtells W90 to write the hopping energies to a file,band_plotdoes the same for the Wannier bands andwannier_plotfor the Wannier orbitals. We usewannier_plot_supercell = 3, 3, 1to prevent Wannier orbitals to be cut off at the edges of the unit cell (extends the unit cell in the xy plane).num_bandsdescribes the number of bands coming from the QEnscfcalculation andnum_wannthe number of Wannier bands to project to.num_iteris the maximum number of iterations W90 is allowed to execute while minimizing the spread of the Wannier orbitals. If the spread is not satisfactory after the given number of iteration, addrestart=wannieriseto the input and W90 will continue where it left off.
W90 disentanglement input
dis_win_min = -12
dis_win_max = 8
dis_froz_min = -7
dis_froz_max = -1.5
dis_num_iter = 500
Description
These arguments set the parameters with which to disentangle the #num_bands to extract #num_wann bands.
The choice of
dis_froz_minanddis_froz_maxis the most important challenge of any W90 calculation. Together they set a window in which only the Wannier bands are present. The window should then be chosen as large as possible for better convergence. In Fig. 16 we clearly see that between energies -7 eV and -1.5 eV only the Dirac bands are present, an easy choice for the disentanglement window.The
dis_minanddis_maxarguments also form a window which should contain all Wannier orbitals and be as small as possible. Fig. 16 indicates that a good choice is the window from -12 eV to 8 eV.dis_num_iteris de maximum amount of disentanglement iterations to perform.
W90 lattice input
begin unit_cell_cart
2.4595121 0.0000000 0.0000000
-1.229756 2.1299999 0.0000000
0.0000000 0.0000000 32.000000
end unit_cell_cart
Begin atoms_cart
ang
C 0.0000000 0.0000000 0.0000000
C 0.0000000 1.4200000 0.0000000
End atoms_cart
Begin projections
C: pz
End projections
Description
These arguments set the lattice and projections.
unit_cell_cartlists the three lattice vectors in Angstrom.atoms_cartlists all atoms and their location in the unit cell.projectionssets the initial guess for the Wannier orbitals as some mix of atomic orbitals. Fig. 16 tells us that a good guess for the Dirac bands are the p\(_z\) orbitals. W90 automatically adds this projection for every atom of the same type. It is also possible to add atomic orbitals to specific locations if we want more control of the exact placement. It is essential that the total amount of projections is identical tonum_wann.
W90 band_plot input
Begin Kpoint_Path
Γ 0.0000000 0.0000000 0.0000000 M 0.5000000 0.0000000 0.0000000
M 0.5000000 0.0000000 0.0000000 K 0.3333333 0.3333333 0.0000000
K 0.3333333 0.3333333 0.0000000 Γ 0.0000000 0.0000000 0.0000000
End Kpoint_Path
Description
The Kpoint_Path is the path in k space W90 follows for the band_plot and has no direct link to the QE calculation.
W90 kpoints input
K points list
mp_grid = 8, 8, 1
Begin kpoints
0.0000000 0.0000000 0.0000000
0.0000000 0.1250000 0.0000000
0.0000000 0.2500000 0.0000000
0.0000000 0.3750000 0.0000000
0.0000000 0.5000000 0.0000000
0.0000000 0.6250000 0.0000000
0.0000000 0.7500000 0.0000000
0.0000000 0.8750000 0.0000000
0.1250000 0.0000000 0.0000000
0.1250000 0.1250000 0.0000000
0.1250000 0.2500000 0.0000000
0.1250000 0.3750000 0.0000000
0.1250000 0.5000000 0.0000000
0.1250000 0.6250000 0.0000000
0.1250000 0.7500000 0.0000000
0.1250000 0.8750000 0.0000000
0.2500000 0.0000000 0.0000000
0.2500000 0.1250000 0.0000000
0.2500000 0.2500000 0.0000000
0.2500000 0.3750000 0.0000000
0.2500000 0.5000000 0.0000000
0.2500000 0.6250000 0.0000000
0.2500000 0.7500000 0.0000000
0.2500000 0.8750000 0.0000000
0.3750000 0.0000000 0.0000000
0.3750000 0.1250000 0.0000000
0.3750000 0.2500000 0.0000000
0.3750000 0.3750000 0.0000000
0.3750000 0.5000000 0.0000000
0.3750000 0.6250000 0.0000000
0.3750000 0.7500000 0.0000000
0.3750000 0.8750000 0.0000000
0.5000000 0.0000000 0.0000000
0.5000000 0.1250000 0.0000000
0.5000000 0.2500000 0.0000000
0.5000000 0.3750000 0.0000000
0.5000000 0.5000000 0.0000000
0.5000000 0.6250000 0.0000000
0.5000000 0.7500000 0.0000000
0.5000000 0.8750000 0.0000000
0.6250000 0.0000000 0.0000000
0.6250000 0.1250000 0.0000000
0.6250000 0.2500000 0.0000000
0.6250000 0.3750000 0.0000000
0.6250000 0.5000000 0.0000000
0.6250000 0.6250000 0.0000000
0.6250000 0.7500000 0.0000000
0.6250000 0.8750000 0.0000000
0.7500000 0.0000000 0.0000000
0.7500000 0.1250000 0.0000000
0.7500000 0.2500000 0.0000000
0.7500000 0.3750000 0.0000000
0.7500000 0.5000000 0.0000000
0.7500000 0.6250000 0.0000000
0.7500000 0.7500000 0.0000000
0.7500000 0.8750000 0.0000000
0.8750000 0.0000000 0.0000000
0.8750000 0.1250000 0.0000000
0.8750000 0.2500000 0.0000000
0.8750000 0.3750000 0.0000000
0.8750000 0.5000000 0.0000000
0.8750000 0.6250000 0.0000000
0.8750000 0.7500000 0.0000000
0.8750000 0.8750000 0.0000000
End kpoints
Description
These k points have to be identical to the k points fo the nscf calculation which is why it is essential to manually define them in the QE input file. mp_grid sets the dimensions of the Monkhorst-Pack grid which has to match the total number of k points.
PW2WAN input¶
The PW2WAN program is part of the QE code which exports the QE wavefunctions for postprocessing by W90.
PW2WAN input
&inputpp
outdir = './out'
prefix = 'graphene'
seedname = 'graphene'
write_unk = .true.
write_mmn = .true.
write_amn = .true.
/
Description
The
outdirandprefixshould match theoutdirandprefixof the QE calculation whereas theseednameshould match the name if the W90 file (minus the extension .win).write_mmnandwrite_amnare crucial inputs for W90 whereaswrite_unkis only necessary forplot_wannier. For a large amount of bands theunkfiles can become massive (couple of GB a piece times the number of k points) so we advise to only add this parameter for small lattices.
Command line interface¶
Before we can run PW2WAN we need to preprocess the W90 input file. So we run the following series of commands:
wannier90.x -pp seedname.win
pw2wannier90.x < prefix.pw2wan.in > prefix.pw2wan.out
wannier90.x seedname.win
W90 offers an mpi version for parellization which has to be compiled seperately. We can parallelize pw2wannier90.x in the same way as pw.x.
W90 output¶
W90 disentanglement output
+---------------------------------------------------------------------+<-- DIS
| Iter Omega_I(i-1) Omega_I(i) Delta (frac.) Time |<-- DIS
+---------------------------------------------------------------------+<-- DIS
1 1.98367895 1.96460709 9.708E-03 0.00 <-- DIS
2 1.97124443 1.95760453 6.968E-03 0.00 <-- DIS
3 1.96331159 1.95227135 5.655E-03 0.00 <-- DIS
258 1.91625911 1.91625911 9.655E-11 0.73 <-- DIS
259 1.91625911 1.91625911 9.135E-11 0.73 <-- DIS
260 1.91625911 1.91625911 8.644E-11 0.73 <-- DIS
<<< Delta < 1.000E-10 over 3 iterations >>>
<<< Disentanglement convergence criteria satisfied >>>
Final Omega_I 1.91625911 (Ang^2)
Description
Omega_I is the gauge invariant part of the spread which minimizes for smooth Wannier bands. Delta is the fractional change of Omega_I:
$\(\Delta = \frac{|\Omega_I(i)-\Omega_I(i-1)|}{\Omega_I(i)}\)$.
When this value drops below the convergence threshold dis_conv_tol the disentanglement process stops.
W90 wannierization output
*------------------------------- WANNIERISE ---------------------------------*
+--------------------------------------------------------------------+<-- CONV
| Iter Delta Spread RMS Gradient Spread (Ang^2) Time |<-- CONV
+--------------------------------------------------------------------+<-- CONV
*----------------------------------------------------------------------------*
Initial State
WF centre and spread 1 ( -0.000000, 0.000000, 0.000000 ) 0.96825063
WF centre and spread 2 ( -0.000000, 1.420000, 0.000000 ) 0.96825057
Sum of centres and spreads ( -0.000000, 1.420000, 0.000000 ) 1.93650120
0 0.194E+01 0.0000000000 1.9365012594 0.73 <-- CONV
O_D= 0.0133595 O_OD= 0.0068826 O_TOT= 1.9365013 <-- SPRD
*----------------------------------------------------------------------------*
*----------------------------------------------------------------------------*
Cycle: 2000
WF centre and spread 1 ( 0.000000, 0.000000, 0.000000 ) 0.96817487
WF centre and spread 2 ( -0.000000, 1.420000, 0.000000 ) 0.96817481
Sum of centres and spreads ( 0.000000, 1.420000, 0.000000 ) 1.93634968
2000 0.000E+00 0.0000000100 1.9363497402 1.85 <-- CONV
O_D= 0.0132080 O_OD= 0.0068826 O_TOT= 1.9363497 <-- SPRD
Delta: O_D= 0.1734723E-16 O_OD= -0.1734723E-17 O_TOT= 0.0000000E+00 <-- DLTA
*----------------------------------------------------------------------------*
Writing checkpoint file graphene.chk... done
Final State
WF centre and spread 1 ( 0.000000, 0.000000, 0.000000 ) 0.96817487
WF centre and spread 2 ( -0.000000, 1.420000, 0.000000 ) 0.96817481
Sum of centres and spreads ( 0.000000, 1.420000, 0.000000 ) 1.93634968
Spreads (Ang^2) Omega I = 1.916259106
================ Omega D = 0.013208013
Omega OD = 0.006882621
Final Spread (Ang^2) Omega Total = 1.936349740
Description
In the wannierization process the gauge dependent spread (Omega_D+Omega_OD, O(D) for (off)-diagonal) is minimized with a steepest descend algorithm ( RMS gradient is the root mean square of the gradient). Delta does not drop below the convergence threshold so W90 completes all num_iter cycles. Our projections work quite well for this system so not much optimization is needed.
Fig. 17 The disentagled wannier bands produces by W90. We compare it to the tight-binding Hamiltonians with NN and NNN hoppings.¶
The NN hopping energy between two Wannier orbitals is approximately 2.9 eV, slightly more than the 2.7 eV of [CDAnjouG+11], which aims to describe the Dirac point more accurately.
Fig. 18 The two Wannier orbitals which closely resemble the p\(_z\) projections.¶